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ABSTRACT 

We present results of a simulation of QCD on a 16 3 x 4 lattice with 2 continuum 
flavors of p4-improved staggered fermion with mass m/T = 0.4. Derivatives of the 
thermodynamic grand potential Q with respect to quark chemical potential \i q up 
to fourth order are calculated, enabling estimates of the pressure, quark number 
density and associated susceptibilities as functions of \x q via Taylor series expan- 
sion. Discretisation effects associated with various staggered fermion formulations 
are discussed in some detail. In addition it is possible to estimate the radius of 
convergence of the expansion as a function of temperature. We also discuss the 
calculation of energy and entropy densities which are defined via mixed derivatives 
of Q with respect to the bare couplings and quark masses. 
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1 Introduction 



Non-perturbative studies of QCD thermodynamics with small but non-zero baryon 
charge density by numerical simulation of lattice gauge theory have recently made 
encouraging progress jT] - jlj. In particular, it has proved possible to trace out 
the pseudo critical line T c (fi q ) separating hadronic and quark-gluon plasma (QGP) 
phases in the (fi q , T) plane out to \i q ~ O(100)MeV j2] - 0j, where the quark chemical 
potential \i q is the appropriate thermodynamic control variable in the description 
of systems with varying particle number using the grand canonical ensemble. In 
addition, the first estimate has been made of the location along this line of the 
critical endpoint expected for iVf = 2 light quark flavors, where the crossover be- 
tween hadron and QGP phases becomes a true first order phase transition |2j. As 
well as being of intrinsic theoretical interest, such studies are directly applicable 
to the regime under current experimental investigation at RHIC, where corrections 
to quantities evaluated at fi q = are both small and calculable. In this respect it 
is worth reminding the reader that in a relativistic heavy ion collision of duration 
~ 10~ 22 s, thermal equilibration is possible only for processes mediated by the strong 
interaction, rather than the full electroweak equilibrium achievable, say, in the core 
of a neutron star. This means that each quark flavor is a conserved charge, and 
conditions at RHIC are thus approximately described by 

jU u = l^d = fi q ; A*i = 2(ju u - /i d ) = 0; /i s = 0, (1) 

with \i q ~ 15MeV [5] when we relate the quark and baryon number chemical poten- 
tials via fiB = 3/x g . In this paper we will present numerical results for the equation 
of state, i.e. pressure p(fi q ,T) and quark number density n q (fi q ,T), obtained from 
a lattice QCD simulation with iVf = 2, which should give a qualitatively correct 
description of RHIC physics, and provide a useful warm-up exercise for the physical 
case of 2+1 flavors with realistic light and strange quark masses. 

Direct simulation using standard Monte Carlo importance sampling is ham- 
pered because the QCD path integral measure det^M, where M(fi q ) is the Eu- 
clidean space fermion kinetic operator, is complex once fi q ^ 0. In the studies 
which have appeared to date, two fundamentally distinct approaches to this prob- 
lem have emerged. In the reweighting method results from simulations at fi q = 
are reweighted on a configuration-by-configuration basis with the correction factor 
[detM(yU g )/detM(0)] 7Vf yielding formally exact estimates for expectation values. In- 
deed, it is found that if reweighting is performed simultaneously in two or more 
parameters, convergence of this method on moderately-sized systems is consider- 
ably enhanced jHj- This method has been used on lattice sizes up to 12 3 x 4 with 
Nf = 2 + 1 to map out the pseudocritical line and estimate the location of the crit- 
ical endpoint at [f q nt ~ 240MeV, T CTtt ~ 160MeV More recently the equation 
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of state in the entire region to the left of the endpoint has been calculated this way 
[7j. However, it remains unclear whether the thermodynamic limit can be reached 
using this technique. 

Analytic approaches use data from regions where direct simulation is possible, 
either by calculating derivatives with respect to fi q (or more properly with respect to 
the dimensionless combination fJ, q /T) to construct a Taylor expansion for quantities 
of interest HU IU E], or more radically by analytically continuing results from 
simulations with imaginary fi q (for which the integration measure remains real) to 
real \i q . The second technique has been used to map T c (fi q ) for QCD with both N{ = 
2 [3] and iVf = 4 ^U] , hi the latter case finding evidence that the line is first order in 
nature. Fortunately, the pseudocritical line found in is in reasonable agreement 
with that found by reweighting; moreover the radius of convergence within which 
analytic continuation from imaginary fi q is valid corresponds to [i q jT < n/3 jllj . 
The leading non-trivial term of quadratic order in the Taylor expansion appears 
to provide a good approximation throughout this region. In general though, while 
analytic approaches have no problem approaching the thermodynamic limit, it is 
not yet clear if and how they can be extended into the region around the critical 
endpoint (but see [12]), and to observables that vary strongly with fi q like eg. the 
pressure or energy density. 

In our previous paper [B] we used a hybrid of the two techniques, by making 
a Taylor series estimate of the reweighting factor [detM(/x 9 )/detM(0)] 7Vf to 0(/i^). 
Since this is considerably cheaper than a calculation of the full determinant, we are 
able to explore a larger 16 3 x 4 system, and also exploit an improved action in both 
gauge and fermion sectors, thus dramatically reducing discretisation artifacts on 
what at T c (fi q = 0) ~ 170MeV is still a coarse lattice. Our results yield a curvature 
of the phase transition line T c (d 2 T c /dfi q )\^ q= o — —0.14(6), consistent with the other 
approaches. Although our simulation employs quark masses m/T = 0.4,0.8, not 
yet realistically light, the results also suggest that any dependence on quark mass is 
weak. In the current paper we extend the Taylor series to the next order O(fi^) but 
this time remain entirely within the analytic framework, using derivatives calculated 
at fi q = to evaluate non-zero density corrections to the pressure p and quark 
number susceptibility Xq = dn q /dfi q , as well as the quark number density n q itself. 
In fact, since the correction Ap can be evaluated at fixed temperature, it turns out 
to be considerably easier to calculate than the equation of state at fi q = [13 EE] . 
Since we now have the first two non-trivial terms in the Taylor expansion, we are 
also able to estimate its radius of convergence as a function of T, and confirm that 
close to T c (fi q = 0) the results of our previous study for the critical line curvature can 
be trusted out to O(100MeV), whereas at higher temperatures a considerably larger 
radius of convergence is likely to be found. Finally we consider mixed derivatives 
with respect to both \x q and the other bare parameters (3 and m, which are required 
to estimate energy e and entropy s densities. Due to the presence of a critical 
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singularity, these latter quantities appear considerably harder to calculate in the 
critical region using this approach. 

Sec. El outlines the formalism used in the calculation and specifies which deriva- 
tives are required. In Sec. El we present a calculation of the cutoff dependence of 
the terms in the expansion for standard, and for Naik- and p4-improved staggered 
lattice fermions, showing that both improvements result in a dramatic reduction 
of discretisation effects. Our numerical results are presented in Sec. 01 and a brief 
discussion in Sec. El Two appendices contain further details on the calculation of 
the required derivatives and the cutoff dependence. 



2 Formulation 

In the grand canonical ensemble pressure is given in terms of the grand partition 
function Z(T, V, fi q ) by 



Note that we have been careful to express both sides of this relation in dimensionless 
quantities. Since the free energy and its derivatives can only be calculated using 
conventional Monte Carlo methods at fi q = 0, we proceed by making a Taylor 
expansion about this point in powers of the dimensionless quantity fi q /T: 



where derivatives are taken at fi q = 0. Note that calculating A(p/T A ) is consid- 
erably easier than p(T,fi q = 0) itself, because whereas InZ must be estimated by 
integrating along a trajectory in the bare parameter plane fSJ [TS], its derivatives 
can be related to observables which are directly simulable at fixed (/3,m), where (3 is 
the gauge coupling parameter and m the bare quark mass. Only even powers appear 
in © because as shown in [3], odd derivatives of the free energy with respect to fi q 
vanish at this point. Note also that we will work throughout with fixed bare mass, 
implying that our computation of A(p/T 4 ) is strictly valid along a line of fixed m/T. 
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For QCD with staggered fermions the partition function may be written 



Z = / VU{debM) N < /A exp{-S g ), (4) 



with U G SU(3) denoting the gauge field variables, S g [U] the link action and 
M[U] f/, q ] the kinetic operator describing a single staggered fermion, equivalent to 
Nf — 4 continuum flavors. On a lattice of size N 3 x iV r with physical lattice spacing 
a, so that T = (aiV T ) _1 , we define a dimensionless lattice chemical potential variable 
fi = ji q a. Equation Q then becomes 

A /p\ liV T 3 2 <9 2 lnZ IN 3 4 <9 4 lnZ 

The derivatives may be expressed as expectation values evaluated at /i = 0: 
<9 2 lnZ _ /JVf d 2 (lndetM) \ / /iV f <9(lndetM) 




a 4 inz 



+3 
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N t <9 4 (ln det M) 
1 d/I 4 



<9 2 (lndetM) 
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NA 2 d 3 (In det M) <9(ln det M) 
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<9/i 3 



<9/i 
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iV f \ 3 d 2 (In det M) / <9(ln det M) 



<9/i 2 



iV f <9(lndet M) 
1 fyx 



-3 



N { <9 2 (ln det M) 
1 d/I 2 



7V f <9(ln det M) 
~4 djl 



1 2 



All expectation values are calculated using the measure Z 1 (fi = 0)T>U (det M[fi = 
])iv f /4 e -s a and in driving © and (J7J) we used the fact that (<9 n (lndetM)/<9/i n ) = 
for n odd. To evaluate these expressions we need the following explicit forms: 



<9(ln det M) 

d/j, 

^ (In det M) 

<9/i 2 
<9 3 (ln det M) 

<9/i 3 



( i dM ' 

;tr r wj 



— M 



+2tr yM~ 
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(8) 
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a 4 (In det M) 



tr M 



r i d 2 M i d 2 M 
-3tr M^-^rM' 1 



4tr ( M" 1 — M" 1 — — 
a/i a/i 3 

12tr M — — M — — M 



dp, 



dp, dp 2 



% , ,3m,, ,9M w ,9M„ ,9M 
-6tr ( M~ — — M~ — — M~ — — M~ — — 
a/i a/i a/i a/i 



The traces can be estimated using the stochastic method reviewed in Since 
d n M/dp n is a local operator, expressions containing p powers of M~ l require p 
operations of matrix inversion on a vector. 

Next we discuss the quark number density n q and its fluctuations. Starting from 
the Maxwell relation 



d 2 tt dN q dp 



dVdp q dV dp q ' 



(12) 



where Q = —Tin Z is the thermodynamic grand potential and N q the net quark 
number, we can write an equation for the quark number density n q analogous to 



T 3 



N 2 d 2 \nZ f IN 2 3 d 4 \nZ 



6 N§ fi dp 4 



+ 



(13) 



It is also possible to interpret derivatives of p with respect to p q in terms of the 
various susceptibilities giving information on number density fluctuations [HI ITT] . 
We define quark number (q), isospin (I) and charge (C) susceptibilities as follows: 



Xq 

rp2 
Xj_ 

J"2 

Xc 

J^2 











d(^/T) 

I ( d_ 

4 \d(fi u /T) 
2 d 



d 



n d 



1 d 



3d(p u /T) 3d(p d /T) 



J"3 ' 

n a - n d 

2n u - n d 
3T 3 ' 



(14) 
(15) 
(16) 



Quark and baryon number susceptibilities are related by xb = dns/dp-Q = 3~ 2 x 9 . 
Any difference between \ q and 4%i is due to correlated fluctuations in the individual 
densities of u and d quarks. With the choice /i u = p d = p q = pa -1 , m u = m d , which 
approximates the physical conditions at RHIC, Xg can then be expanded in terms 
of quantities already used in the calculation of p and n q . 



Xq_ 
J^2 



1 dn q N T d 2 \\\Z_ d 2 ( Xq /T 2 



M 9 =0 



T 2 dp q iV 3 dp 2 



d(N/T) 2 



d 4 lnZ 



M 9 =0 



N T N% dp 4 



(17) 
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whereas the expansion of Xi 1S determined by the following expectation values: 



Xi 

d\ X i/T 2 



N T /2<9 2 (lndetM) 
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2 \ 2 <9 3 (In det M) d(ln det M) 



dfi 3 



dfi 



2\ 3 <9 2 (lndetM) /d(lndetM) 



2 <9 2 (In det M) 
4 ^u 2 



2 <9(ln det M) 
4 



dfi 

2 a 2 (In det M) 
4 iV 



(19) 



where we have explicitly set Nf = 2. Charge fluctuations are then given by the 
relation 



X° 1 X g Xi 1 ( d{n n /T*) d(n d /T 3 ) \ 
T 2 36T 2 T 2 6 V <9(/i u /T) <9(/i d /T) / 

where the third term vanishes for fi u = /id, Tt u = m^. 



(20) 



Finally we discuss the energy density e, most conveniently extracted using the 
conformal anomaly relation 



3p 



VT 3 



0(3 d In Z dm d\nZ 

da dp da dm 



(21) 



where (3 and m are the bare coupling and quark mass respectively. In fact, for fi ^ 
the derivation of this expression needs careful discussion. Start from the defining 
relation 



Q = E - TS - fi q Nq = -pV = -TlnZ, 



(22) 



where S is entropy. For a Euclidean action S = S(/3, m, /i) defined on an isotropic 
lattice of spacing a we have the identity 



dS Tr dS dS 

= 3V T — . 

dV dT 



da 

It follows that 
dQ 



V 



dV 



VT 



dS 
dV 



dtt 
T — 
dT 



. dS 
<9T 



(23) 



-pV 

= -TS = Q-E + f.i q N q 



(24) 
(25) 
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implying 



e-3p-fi q n q 



T 
V 



dS_ 

da 



T 
V 



dp din Z dm dlnZ d/jdlnZ 
da dp da dm da dji 



(26) 



where we have allowed for the dependence of the lattice action on all bare parameters. 
Since however /i = [i q a, and a parameter multiplying a conserved charge experiences 
no renormalisation, the third terms on each side cancel leaving the relation (J21)) . 

Taylor expansion of (J21j) about fi = leads to the expression 



3p 
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dpN* 
l da~Nl 
dm N 3 
daPP 



,d 3 InZ 
dpdfi 2 ' 
d 3 \nZ 



24^ 



<9 5 \nZ 
dpd^ A 
d 5 \nZ 



+ 



1 2 d 6 \nZ 1 
2^ dmdp 2 + 24^ dmd/j, 4 



(27) 



The beta function a(dp/da) may be estimated by measurements of observables at 
(T,fi q ) = (0,0); the factor a(dm/da) is poorly constrained by current lattice data 
but vanishes in the chiral limit, so is frequently neglected. In order to assess the 
magnitude of the resulting error, it is nonetheless useful to calculate all the derivative 
terms. They may be estimated using the formulae 



d(O) 
dp 

d(Q) 
dm 



o 

dO 

dm 



dSg 

' dp 



(O) 



OSr, 



O 



dp 

Nf <9(lndetM) 
4 dm 



(O) 



Nf <9(lndet M) 
4 dm 



(28) 
(29) 



The derivative dS g /dp is, of course, simply the combination of plaquettes comprising 
the gauge action itself, and derivatives with respect to m can be evaluated using 



<9 n+1 (lndetM) ^(trM" 



dmdfi r 



d/j, n 



(30) 



The implementation of the second square bracket in (|27|) in terms of lattice operators 
is straightforward but unwieldy; for reference the non-vanishing terms are listed in 
Appendix 1X1 



3 Analyzing the cut-off dependence 

In this section we discuss the influence of a non-zero chemical potential n q on the 
cut-off effects present in calculations of bulk thermodynamic observables on a lattice 
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with finite temporal extent N T . For \i q = this issue has been discussed extensively 
for both gluonic and fermionic sectors of QCD. In particular, it has been shown that 
the use of improved actions is mandatory if one wants to ensure that discretisation 
errors in the calculation of quantities like the pressure p or energy density e are 
below the 10% level on moderately sized lattices N T <(8 — 10) jHj. We now want to 
extend these considerations to the 7^ 0, which affects the quark sector only. 

Following [T3| we will concentrate on an evaluation of the pressure. As we will be 
evaluating thermodynamic quantities using a Taylor expansion in fi q /T we want to 
understand the cut-off dependence of p(fi q ) and its expansion coefficients in terms 
of n q /T. 

In the limit of high temperature or density, due to asymptotic freedom ther- 
modynamic observables like p or e are expected to approach their free gas, i.e. 
Stefan-Boltzmann (SB) values. In this limit cut-off effects become most significant 
as the relevant momenta of partons contributing to the thermodynamics are 0(T) 
and thus of similar magnitude to the UV cut-off a -1 . Short distance properties thus 
dominate ideal gas behaviour and cut-off effects are controlled by the lattice spacing 
expressed in units of the temperature, Ta = 1/N T . 



by 



In the continuum the pressure of an ideal gas of quarks and anti-quarks is given 



p_ 

y4 



Nf 



POO 

/ dkk 2 \n [(1 + zexp{-e(k)/T})(l + z~ l exp{-e(k)/T})] (31) 
Jo 



2tt 2 T 3 

with the fugacity z = exp{fi q /T} and the relativistic single particle energies e(k) = 
\/k 2 + m 2 . For massless quarks one finds from an evaluation of the integral the 
pressure as a finite polynomial in fi q /T: 



P_ 



7N { 7i 2 Nf /M 2 Nj_ //M 



60 +iH?J +sH?J ■ < 32 > 



For m non-zero the pressure is a series in the fugacity: 

2 AT °° 

^£(-l) m ^(WT)(/ + ^), (33) 

where K 2 is a Bessel function. Of course, Eq. (}3*3*j) can also be reorganised as a 
power series in \i q jT . 

It is well known that the straightforward lattice representation of the QCD par- 
tition function in terms of the standard Wilson gauge and staggered fermion actions 
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leads to a systematic 0(a 2 ) cut-off dependence of physical observables. In the infi- 
nite temperature limit this gives rise to 0((aT) 2 = 1/N 2 ) deviations of the pressure 
from the SB value (j32j) ; 



p_ 

rp4 



J)_ 

rp4 



(34) 



Using improved discretisation schemes it is possible to ensure that cut-off effects 
only start to contribute at 0(N~ 4 ) [T3], or to considerably reduce the magnitude of 
the coefficient d relative to the standard discretisation scheme for staggered fermions 

For fi q = the pressure of free staggered fermions on lattices with infinite spatial 
volume (N a = oo) but finite temporal extent N T is given by 



p_ 

rp4 



8 f T (2n)3j 



1 f 27T 

-— J dp 4 ln(u 2 (p) + 4f 2 (p)) 



2- 



d 3 p 



iV-^ln (o; 2 (p) + 4/ 4 2 (p)) 



(35) 



In the first term the sum ^ runs over all discrete Matsubara modes, i.e. p^ G 
{(2n + l)ir/N T \n = 0, . . . , N T — 1}, whereas in the second term we have an integral 
over p 4 which gives the vacuum contribution. For quarks of mass m the function 
u 2 (p) is given by u 2 (p) = 4^^ =1 / 2 (p) + N~ 2 (m/T) 2 . Here we have introduced 
functions to describe the momentum dependence of the propagator for the 

standard, Naik jTSj and p4 staggered fermion actions 



Up) 
Up) 
Up) 



- sin p u 

9 . 1 . „ 

— sm Pll - — sm3p M 



(standard staggered action) 
(Naik action) 



| sin p^ + — sin p^ ^ cos 2p u (p4 action) . 



(36) 
(37) 
(38) 



The introduction of a non-zero chemical potential is easily achieved by substituting 
every temporal momentum p 4 by p^ — %\i = p^ — zA r ~ 1 (/i g /T). The integrals in (J33|) 
have been evaluated numerically for different N T . Results for different values of [i q jT 
and m/T are shown for the different fermion actions in Fig. ^ For the standard 
action cutoff effects remain > 10% out to N T pa 16, whereas both improved actions 
are hard to distinguish from the continuum result by N T = 10. We note that lines 
for different [i q /T values but the same quark mass fall almost on top of each other. 
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2 4 6 8 10 12 14 16 18 20 2 4 6 8 10 12 14 16 18 20 



(a) (b) 
Figure 1: The pressure calculated on lattices with temporal extent N T in units of 
the continuum ideal Fermi gas value, (a) shows results for the standard, Naik and 
p4 actions at (yU g /T, m/T) = (0, 0), (0, 1), (1, 0) and (1, 1); (b) the coefficients 
Co, C2, C 4 of the fiq/T expansion of p{m/T = 0) divided by the corresponding SB 
constant as a function of N T . 



Cutoff effects are thus almost independent of fi q . The effect of \x q 7^ on the cutoff 
dependence of the pressure is even smaller than the effect of quark mass m^0. 

As can be seen from (J32*j) for moderate values of fi q /T the /i-dependence of the 
continuum ideal gas pressure is dominated by the leading 0((/i q /T) 2 ) contribution. 
In order to control the cut-off dependence of the various expansion terms we have 
expanded the integrand of (J55)l up to order 0((/i g /T) 6 ). For the standard action the 
series starts with 



In ( u 2 (p) + sin 2 ( p A - 



N T T 
2i cosp4 sinp 4 r\i q 
DN T V ~f 
— 1 + 4 D cos 2p4 + cos 4p 4 / ji q \ '■ 

i (— 1 + 4 D 2 + 6 D cos 2p 4 + cos 4p 4 ) sin 2p 4 f^ q \ 3 
6 D 3 N* V ~f 



+ O 



(39) 
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Here we use the shorthand notation D = 4 Ylu=i f^(p)- The remaining orders as well 
as the series for Naik and p4 actions are given in Appendix [Bj A common feature 
of these expansions is that the odd terms are pure imaginary and the integral and 
sum over p 4 of those terms vanish due to a factor sin (77.^4) which always appears. 
To be more precise, this factor always forms the pattern sin(np4) cos(mp4) which 
can be shown to vanish, either after summation over the discrete set of values, 
or integration from to 2tt, for n, m G N. Performing the momentum integration 
and the summation over Matsubara modes we obtain the coefficients of the ix q /T 
expansion of the pressure; 



P_ 



We checked numerically that with increasing N T the coefficients Co, C2 and C4 do 
indeed converge to their corresponding SB values, 

lim Co = - — ; lim C2 = -; lim C4 = — -. (41) 

N T ^oo 60 N T -*oo 2 N T ->oo An 2 

Fig. shows Co, C2 and C4 for the standard, Naik and p4 actions with massless 
quarks, normalized to the corresponding SB value. We see here again that the cutoff 
dependence of the pressure at fi 7^ is qualitatively the same as at /x = 0. 

For massless quarks the coefficient Cq should vanish with increasing N T , as 
checked in Fig. |2K It is expected that this term will approach zero like N~ 2 in 
the large N T limit. In order to define the numerical factor, we plot CqN 2 over N~ 2 . 
A fit yields C$ ~ — 0.015 N~ 2 for the standard action. This is at least an order of 
magnitude larger than for the p4 improved action, for which the dominant cut-off 
dependence seems to be 0(N~ 4 ) as for the Naik action. 

In the case of massive quarks the expansion (|4()jl no longer terminates at O(fig). 
After expanding ()31|) in terms of n q /T and performing a numerical integration we 
find for the expansion coefficients Ci(m/T) up to % = 6 the values given in Tabled 
The mass value m/T — 0.4 is the value we use in our numerical calculations, cor- 
responding to N T = 4 and am = 0.1. We note that the coefficient Cq no longer 
vanishes. As shown in Fig. |2b ; for iV T finite there are large deviations from the 
continuum value. Even at N T = 4, however, the absolute value of this coefficient is 
still a factor of about 10~ 4 smaller than the leading term Co. The deviations thus 
do not show up in the calculation of the complete expression for the pressure shown 
in Figure HJi. These terms, however, become more important in higher derivatives 
of the pressure such as the quark number susceptibility Xg- I n summary, for a gas 
of free quarks we find that the fi q /T expansion up to 0((/i g /T) 4 ) is sufficient for 
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(a) (b) 
Figure 2: (a) The coefficients Cq in the massless case, multiplied with iV 2 as a 
function of N~ 2 (results for standard staggered fermions are divided by 10) and (b) 
the ratio Cq(N t )/Cq(oo) at m/T = 0.4 for the standard, Naik and p4 actions. 



\x q jT < 1 and m/T < 1. In the continuum the deviation from the full expression 
over this range is smaller than 0.01%. On the lattice, however, cut-off effects lead 
to deviations of approximately 10% on coarse (N T = 4) lattices. 



4 Numerical Results 

We applied the formalism outlined in Sec. El to numerical simulations of QCD with 
Nf — 2 quark flavors on a 16 3 x 4 lattice, using both Symanzik improved gauge 





m/T 


= 0.4 


m/T 


= 1.0 


i 


Ci{m/T) 


Ci(m/T)/Ci(0) 


Ci{m/T) 


Ci(m/T)/Ci(0) 





1.113632 


0.967 


9.528163- 10" 1 


0.827 


2 


4.880455 • 10" 1 


0.976 


4.313914- 10" 1 


0.863 


4 


2.531101 • 10~ 2 


0.999 


2.471397- 10~ 2 


0.976 


6 


1.877659 • 10" 6 




5.036816- 10" 5 





Table 1: Continuum values for the coefficients C\ of the /x 9 /T expansion of the 
pressure of a massive gas of quarks for the mass values m/T = 0.4 and m/T = 1.0. 
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and p4-improved staggered fermion actions. The simulation method is exactly as 
presented in |3j a . The bare quark mass was ma = 0.1 for which the pseudocritical 
point for zero chemical potential is estimated to be f3 c ~ 3.649(2). In order to cover 
a range of temperatures on either side of the critical point we examined 16 values 
in the range (3 G [3.52,4.0]. The simulation employed a hybrid molecular dynamic 
'R'-algorithm with discrete timestep St = 0.025, and measurements were performed 
on equilibrated configurations separated by r = 5. In general for each (3 500 to 800 
configurations were analysed, with 1000 used in the critical region (3 G [3.52,3.66]. 
On each configuration 50 stochastic noise vectors were used to estimate the required 
fermionic operators. For each noise vector, 7 matrix inversions are required to 
estimate the required operators (jHHUJ) and (j49H52|) . 

Following the procedure used for the equation of state at n — [ISj, we translate 
to physical units using the following scaling ansatz [TSj : 



where a = R{(3)/R{(3), R being the two-loop perturbative scaling function appropri- 
ate for two light flavors. Using string tension data at (T, //) = (0, 0) the best values 
for the fit parameters corresponding to a reference (3 = 3.70 are g<i = 0.669(208), 
<74 = —0.0822(1088) We find that our simulations span a temperature range 

T/T G [0.76, 1.98], where T is the critical temperature at \x q = 0. 

In Fig. El we show the first two coefficients, (a) c 2 and (b) C4, of the Taylor 
expansion of A(p/T 4 ) introduced in (J3J) as functions of T/Tq. Also shown are the 
corresponding SB limits: (a) N{C2(N T ) and (b) N{C4(N T ), where the coefficients Ci 
are defined in (|40[). with values relevant for both the lattice used (N T = 4) and 
the continuum limit (N T = 00) plotted. Both c 2 and C4 vary sharply in the critical 
region, but except in the immediate vicinity of the transition the quadratic term 
dominates the quartic. This is consistent with the results of [7j where data at 
varying fi obtained by reweighting was found lie on an almost universal curve when 
plotted as a fraction of the SB prediction. The asymptotic value of C4 appears to be 
approached from above. 

A notable feature is that in the high-T limit our data lies closer to the continuum 
SB prediction rather than their values Ci(N T = 4) corrected for lattice artifacts, C2 
assuming 80% of the continuum value for T/T = 2 whereas c 4 is almost coinci- 
dent with its continuum value. By contrast recent calculations with unimproved 
staggered fermions [7| find that the high-T limit of the data lie close to the 
lattice-corrected SB value. This situation can be modelled by making the coefficient 

a Thc coefficient cf of the knight's move hopping term was incorrectly reported to be 1/96 in 
[2; its correct value is 1/48. 



l + 92a 2 (f3)+g^m 

1+92+94 



(42) 



11 



0.2 



1 1 1 ' 


1 1 1 1 1 1 1 




SB limit 




_— — e— -© 9 ? 




SB (AT=4)" 


&T ,1,1 


1,1,1, 



0.25 



0.8 1 1.2 1.4 1.6 1.8 2 

r/r 

(a) 




0.15 - 



0.05 



77T„ 



(b) 



Figure 3: Coefficients of (a) (p q /T) 2 and (b) (/i 9 /T) 4 in the Taylor expansion of 
A(p/T 4 ) as functions of T/Tq. 



d of the 0{N~ 2 ) correction appearing in (|34|) temperature dependent. In thermody- 
namic calculations performed with pure unimproved SU(3) lattice gauge theory [T9] . 
where extrapolations to the continuum limit are currently practicable, it is found 
that d(T) ~ 0.5d(T = oo) for T ~ 3T , becoming even smaller closer to T . The 
behaviour of C2 and C4 we have observed using p4 fermions is broadly consistent with 
this behaviour. 

In Fig. |U we plot A(p/T 4 ) defined in (J3J) as a function of (fi g /T) 2 for various 
temperatures. In most cases C4 <C C2 and the relation thus almost linear. The 
strongest departures from linearity are for T ~ T , but even here the quadratic term 
is dominant for (p, q /T) 2 <0A, corresponding to /i g <100MeV. Given enough terms of 
the Taylor expansion in fi q /T, one could determine its radius of convergence p via b 



p = lim p r , 



lim 

n— >oo 



Cn+2 



(43) 



Data from the pressure at p q = [T^] and the current study enable us to plot the first 
two estimates po and P2 on the (p q , T) plane along with the estimated pseudocritical 
line T c (p q ) found in j3] in Fig. El Also shown are the corresponding values from the 

b The argument of @] that p < f due to the presence of a phase transition as imaginary chemical 
potential is increased beyond this value [201 does not hold for calculations with fx q real; in this case 
the pressure po corresponding to the unit element of the Z(3) sector is always the maximum and 
hence dominates the partition function in the thermodynamic limit - hence the issue concerns the 
analytic properties within this physical unit sector. 
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Figure 4: A(p/T A ) as a function of (n q /T) 2 for various temperatures, increasing 
upwards from the lowest curve with T/Tq = 0.812 to the highest with T/Tq = 1.980. 
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SB limit (132)1 . For T > T c one finds that p n increases markedly as n increases from 
to 2; if the SB limit is a good predictor for the QGP phase we might expect c 6 to 
be very small, and the next estimate p4 correspondingly very large in this regime. 
Close to the transition line, however, the thermodynamic singularities appear to 
restrict p ~ 0(1); this in turn gives an approximate lower bound for the position of 
the critical endpoint. From the figure we deduce p q Tlt >(l — 1.2)To, not inconsistent 
with the result of The new results at 0(/i 4 ) are important because they justify in 
retrospect our neglect of fourth order reweighting factors in our earlier calculation 
of the critical line T c (p q ) [3]. Indeed, simulations with imaginary p suggest that 
neglect of these terms in the analytic continuation to physical \i q is justified for 
^<170MeV @j. 

Fig. El plots the dimensionless correction A(p/T 4 ) to the equation of state as 
a function of both p q /T and \i q . In the latter case the correction rises steeply 
across the transition and peaks for T ~ 1-1T , before rapidly approaching a form 
A(p/T 4 ) = aT~ 2 characteristic of the SB limit, with the coefficient a having 82% 
of the continuum SB value. Comparison with the equation of state results at \i q = 
from Ref. (T^j suggest that the correction will give a significant correction to the 
pressure for 0.9<T/T <1.3, fi q /T >0.5, but will decrease in importance as T rises 
further. The curves of Fig. \§jp are in good qualitative agreement with those of 
Refs. j?J E] , although we consider any quantitative agreement to be somewhat acci- 
dental as the numerical data obtained in |5| with unimproved actions have large 
discretisation errors which have been corrected for by renormalising the raw data 
with the known discretisation errors in the infinite temperature limit. Experience 
gained in calculations of thermodynamic quantities in the pure SU(3) gauge theory 
suggests that in the temperature range of a few times T this procedure overestimates 
the importance of cut-off effects by a factor two or so ill)] . 

Fig. shows the quark number density n q evaluated using Eqn. (|T3|) . As p q 
increases, n q rises steeply as the QGP phase is entered; for reference, if the quark 
number density in nuclear matter is denoted n q , then the ratio Uq/T^ rs 0.75. Our 
results are numerically very similar to those obtained using exact reweighting in [2] , 
where a mass ma ~ 0.1 for the light quark flavors was used. Note that a significant 
quark mass dependence for n q was observed in [3] , and indeed is present even in the 
SB limit as described in Section El however analysis of the SB limit suggests that 
the difference between the chiral limit and m/T = 0.4 is about 4%. In Fig. we 
show the result of eliminating p q in favour of n q via 
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Figure 7: (a) n q /T 3 as a function of T/Tq for various fi q /To, and (b) the "true" 
equation of state A(p/T 4 ) vs. (n q /T 3 ) 2 for various temperatures. The continuum 
SB forms \(n 2 /N^injT 3 )^ 3 (low T) and (2N i )- 1 (n q /T 3 ) 2 (high T) are also 
shown as functions of T/Tq. 
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The relation (|44j) approximates the "true" equation of state in terms of physically 
measurable quantities; we have plotted the resulting A(p/T 4 ) against (n q /T 3 ) 2 up 
to the point where the ratio of the magnitude of the second term of ()44j) to that 
of the first is 40%: the point n q /T 3 = y^^cfT^Q where the ratio is 50% marks a 
mechanical instability dp/dn q = 0, which is an artifact due to the truncation of 
the series. Stability of the equilibrium state under local fluctuations 5n q requires 
dp/dn q > 0, an example of Le Chatelier's principle. As T/T increases through 
unity, the equation of state changes from a form resembling the low-T SB limit 
p oc n q ^ 3 to the stiffer p oc n 2 q characteristic of the high-T SB limit. Interestingly 
enough, to the order we have calculated the instability artifact sets in at [i q jT ~ 1.4 
for T large, but at fi q /T ~ 0.4 for T « T , thus providing an independent, and 
more stringent, limit to the physical validity of our approach, and reflecting the 
importance of contributions from higher orders in the Taylor expansion close to 

Next, in Fig. |S]we plot the expansion coefficients corresponding to the various 
susceptibilities defined in fll4Hlfjj) . For T<T Q there is a significant difference between 
Xqi^q — 0) and 4xi(/i g = 0), implying anti-correlated fluctuations of n n and which 
rapidly decrease in magnitude above T and vanish as T approaches the infinite 
temperature SB limit c . In the same limit the charge susceptibility xc approaches 

c There has recently been a discussion whether the difference (4%/ — Xq)/T 2 is exactly zero in 
the high temperature phase, as suggested by some lattice calculations or just small but non- 
zero, as suggested by HTL-resummed perturbation theory We find that the difference stays 
non-zero but decreases by one order of magnitude between T ~ To and T ~ 1.5Tb. At T ~ 1.36Xb 
we find a value of 0.0066(28) for this difference calculated in 2-fiavor QCD which clearly disagrees 
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T/T n 



Figure 9: Xq/T 2 as a function of T/T for various ji q /T. 



the value j$Xq- The critical singularity in 4xi and xc is weaker than that of Xq, which 
can be traced back to the differing coefficients of ((<9 2 lndetM/<9/i 2 ) 2 ), the dominant 
term in the vicinity of T c , in the definitions (J7J) and ()19|1. The dimensionless quantity 
Txd s , where s = (e + p — /i q n q )/T is the entropy density, can be related to event- 
by-event fluctuations in charged particle multiplicities in RHIC collisions, and has 
been proposed as a signal for QGP formation [22] • Event-by-event fluctuations in 
baryon number have also recently been discussed in 



In Fig.|21the relation (|17[) and data of Fig.|3]have been used to plot the dimension- 
less quark number susceptibility Xq/T 2 as a function of T/T for various /i q /T. The 
peak which develops in Xq as \i q increases is a sign that fluctuations in the baryon 
density are growing as the critical endpoint in the (/i, T) plane is approached. Phys- 
ically, this shows that at the critical point, as well as strong fluctuations in the (ijjijj) 
bilinear expected at a chiral phase transition there are also fluctuations in (ifj'joijj) 



with the quenched result (2±4) • 10 -6 presented in as well as the recent 2-flavor results of this 
group [3J. Our results are, however, in agreement with the findings of Ref. At T ~ 2Tq the 
numerical value of this difference drops below our current error level of about 2 • 10~ 3 . In the high 
temperature limit this error is thus not yet small enough to discuss numerical effects at the level 
of 10 -4 as suggested in the discussion presented in |23| . 
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Figure 10: Derivatives (a) (2VT 3 )~ 1 d 3 In Z/d/3d(/i q /T) 2 (circles) and 
-(2VT 3 )- 1 d 3 \nZ/dmd(fi q /T) 2 (diamonds); (b) (2AVT 3 )~ 1 d 5 \nZ/d(3d(fi q /T) 4 
(circles) and -(2WT 3 )~ 1 d 5 In Z/dmd(fi q /T) 4 (diamonds). 



since Lorentz symmetry is explicitly broken by the background baryon charge den- 
sity. For quantities such as n q and Xq defined as higher derivatives of the free energy 
with respect to fi q , the relative importance of the higher order terms in the Taylor 
series expansion is increased; for example, at T ~ T and n q /T = 1 the quadratic 
contribution to x<?(/A?) is about 3 times that of the leading order term. For this rea- 
son we do not expect the data of Fig. El to be quantitatively accurate in the critical 
region. Note, however, that at each temperature the expansions for p, n q and \ q all 
have the same radius of convergence. 

Finally we turn to a discussion of the derivatives necessary for calculating the 
response of the energy density e to increasing fi q . The Taylor expansion of the energy 
density involves derivatives of the expansion coefficients c p (T) used to calculate the 
pressure, 

*{^)=p r m(^y • («) 

with c' p (T) = T(dc p (T)/dT)\ flq=0 . It is apparent from the temperature dependence 
of the expansion coefficients C2(T) and c±(T) shown in Fig. El that the coefficients 
c' p (T) can become large in the vicinity of T . On the other hand that figure also 
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shows that c' p (T) will be small, i.e. close to zero, at high temperature as expected 
in the ideal gas limit. A comparison with (J27)l shows that the numerical evaluation 
of c' p (T) requires the knowledge of lattice beta-functions and a calculation of mixed 
derivatives of In Z with respect to \i as well as (3 and m. In Fig. we plot these 
derivative terms; the signals in this case are much noisier than for d n In Z/dfi n , 
although we have been able to check that the numerical values for d 3 In Z/dj3dfi 2 
are consistent with the slope of the curve in Fig. It is clear firstly that with 
the exception of d 3 In Z/dmdfi 2 the signal only differs significantly from zero in 
the immediate neighbourhood of the transition, and secondly that derivatives with 
respect to m are strongly anti-correlated with those with respect to (3. The latter 
suggests it might be possible to learn something from (J27|) about the shape of the 
A((e — 3p)/T 4 ) curve as a function of T/T away from the chiral limit even in the 
absence of quantitative information about adm/da. 

Consider however ignoring mass derivatives and focusing on those performed with 
respect to coupling. In this case all derivatives are consistent with zero for T>1.2 To; 
i.e. the difference A ((e — 3p)/T 4 ) is to a good approximation independent of \x q for 
these high temperatures. This observation is consistent with the results obtained 
by using exact reweighting in [7j. Consider now temperatures close to T . The beta 
function at the critical f3 c has the value a~ 1 da/df3 = —2.08(43) [T^j; substituting the 
derivatives from Fig. El i n ()27|) we find at To 

A (^f?) «(4.8±1.2)x (^) 2 -(5±4)x (^)V-- (46) 

Taking the central values of the coefficients in this expansion one may conclude that 
the ratio c' 2 /c' A is comparable with C2/C4. At present the large error on the coefficient 
of the (/iq/T) 4 term, however, does not allow a firm conclusion on the convergence 
radius of the expansion of e — 3p. We also note that the coefficient c f 4 will change sign 
for T ~ To. This suggests that large cancellations can occur for fi q /T ~ 0(1) and 
indicates that higher order terms are needed to determine this difference reliably. In 
any event, it would appear that extending our current analysis to determine energy 
and entropy densities (e, s) in the critical region will be far more demanding. 

5 Summary 

We have presented the first Monte Carlo calculation of the QCD equation of state at 
non-zero quark chemical potential within the analytic framework; no reweighting has 
been performed. As in our previous work, we have exploited the relative simplicity 
of the method to explore larger physical volumes than those used in comparable 
studies Ej. In addition, the compatibility of our method with the use of an 
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improved lattice fermion action has meant that our results suffer from relatively 
mild discretisation artifacts, our data for the pressure correction Ap(fi q ) achieving 
80% of the Stefan-Boltzmann value by T ~ 2Tq. 

Our results for Ap and its //-derivatives n q and the various susceptibilities \% are 
in good qualitative agreement with those of [7J E] . Since higher derivatives suffer 
from larger discretisation artifacts, and are inherently noisier when estimated by 
Monte Carlo simulation, the results for, say, Xq are l ess quantitatively reliable than 
those for Ap; nonetheless the singularity developing in Xg as [x q is increased, seen in 
Fig. El is evidence for the presence of a critical endpoint in the T) plane, and 
for the importance of quark number fluctuations in its vicinity. 

The calculation of fourth order derivatives has enabled us for the first time to 
make an estimate for the limitations of our method, both analytically through the 
radius of convergence of the Taylor expansion in fi q /T, and physically via the re- 
quirement of consistency with Le Chatelier's principle. For both criteria the most 
stringent bounds are unsurprisingly in the vicinity of the critical line, where conver- 
gence of the series limits us to fi q /T<l, and mechanical stability of the equilibrium 
state to fi q /T<0.5. Of course, the picture should change with the inclusion of still 
higher derivatives since on physical grounds we expect stability of the equilibrium 
state everywhere within the domain of convergence. We are currently investigating 
the feasibility of including the relevant 0(/i 6 ) terms in our calculation. 

Other quantities of phenomenological importance such as the energy e and en- 
tropy s densities, which require mixed derivatives with respect to the other bare 
parameters (3 and m, appear more difficult to calculate with quantitative accuracy 
with this approach. It remains an open question whether the radius of convergence 
for these quantities is the same as that for quantities defined by series in d n Q/dfi n . 

Finally, it is necessary to stress the importance of refining the current calculation, 
firstly by simulating systems with N T > 6 so that a reliable extrapolation to the 
continuum can be performed, and secondly by repeating it with a realistic spectrum 
of 2+1 fermion flavors. 
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A Derivatives needed to calculate energy density 



Here we present the non- vanishing terms in the expressions for d n+ In Z/dmdji n : 



d 3 \nZ 
dji 2 dm 
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dji^dm 
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+ 12 
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+4 
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Nf dhrM- 
1 ~djf 
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'JV f y 9(lndet M) <9trM _1 
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As explained above, all terms involving the expectation value of an odd number of 
derivations with respect to n have been set to zero. Evaluation of eqns. (J47)) and 
(J48J) requires the following expressions for the derivatives of trM -1 : 



dtrM" 1 / id-M" A 

— = " tr { M w M ) (49) 



d/j, \ d/j, J 

a 9 = -tr M" 1 — — M" 1 + 2tr M — — M — — M 
on \ on J \ on on j 

o 3 = "tr M' 1 — -M" 1 + 3tr M — — ^-M — — M , 

Of! 6 \ Ofi J \ Ofi OfJ, J 

+3tr M" 1 — - M~ l — -M~ l - 6tr M ^—M ^—M -—M 

\ on on j \ on on on 

dHrM- 1 ( ,<9 4 M A A ( l d 3 M 1 dM A 
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(50) 
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-12tr M" — -^-M" — — M" — — M~ 
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-12tr M" 1 — - M" 1 — — M" 1 — - M" 1 
\ a/i dji 2 o/j, 

-12tr M — — M — — M — 
\ ctyx o\i o 



2 M 



■M 



\ d/j, dji <9/i 2 J 

+24tr M" 1 ^— M" 1 ^— M" 1 ^— M" 1 ^— M" 1 
\ a/i a/i a/i a/i / 

B The pressure of free staggered fermions 

Expanding the quantity p/T 4 as discussed in section |3]one finds 
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(53) 



Here only the even expansion coefficients give non-vanishing contributions. Intro- 
ducing the abbreviation, 



(54) 



with fn(p) as given in (jHEJ), the even expansion coefficients for the standard action 
are given by: 
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96L> 4 

+ (12 - 56L> 2 ) cos(4p 4 ) - 24Dcos(6p 4 ) - 3 cos(8p 4 )) 
(150 - 180D 2 + 32D 4 - 8,0(45 - 60D 2 

+ 16L> 4 ) cos(2p 4 ) + (-225 + 960L> 2 - 992£> 4 ) cos(4p 4 ) 
+540L>cos(6p 4 ) - 1440D 3 cos(6p 4 ) + 90cos(8p 4 ) 
-78(LD 2 cos(8p 4 ) - 180-Dcos(10p 4 ) - 15 cos(12p 4 )) , 



(55) 
(56) 

(57) 



(58 
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For the Naik action we introduce an additional function, 



g 4 (p) = -%- 



The even expansion coefficients can then be written as: 



9 3 
= cos (P4) + 7^ cos(3p 4 ), (59) 



c 2 = 



c 4 = 



c 6 = 



Co = HD) (60) 
g=|( - QDfj(p) + 6Dgj(p) - 48f 4 (p)gl(p) 

+Df A (p)sm(3p A )) (61) 

3-^4 (48 ( - 768m P )gt(p) + £ 3 (/ 4 » - g 2 M) 
-lV2Dfl(p)gi(p)(ft(p) - glip)) + D 2 ( - 6/ 4 (p) 
+44/ 2 (p)< ?4 2 (p) - 6gt(p))) - 2AD 2 {D - 8/|(p))^(p) cos(3p 4 ) 
—32Df A (p) [D 2 - 3Df 2 (p) + 9Z^(p) - 48/ 2 (p)< ?4 ») sin(3p 4 ) 

+D 2 (D-8/ 2 (p))sin 2 (3p 4 )) (62) 

720D 2 g 4 (p)(D* + 768ft(p)gj(p) + D 2 ( - 2Qf 2 (p) 

+6g 2 4 (p)) + 96D(ft(p) - 2fl{p)gl(p))) cos(3p 4 ) 
-A5D 4 (D - 8/£(p)) cos 2 (3p 4 ) - 96Df 4 (p) (8D 4 + 46080/ 4 4 (p)^ 4 4 (p) 
— 7hD (fl(p) - 3gl(p)) + 60L> 2 (3/ 4 4 (p) - WflMip) + ^gjip)) 
+2880D(3ft(p)g 2 M - 5/|(p)^(p))) sin(3p 4 ) + 15 J D 2 (5 J D 3 
+4608/ 4 (p)s 2 (p) + D 2 ( - 76/ 2 (p) + 36(7 2 (p)) + 192D(ft(p) 
-Qfl(p)gl(p))) sin 2 (3p 4 ) + 10D 3 U(p)(3D - 16/ 2 (p)) sin 3 (3p 4 ) 
+72(4(245760/ 4 (p) 6 <? 4 6 (p) + D 5 (f 2 (p) - g 2 {p)) 
+921Q0Dft(p)gt(p)(fl(p) - glip)) - 2L> 4 (l5/ 4 (p) - 94/ 4 2 (p)« ?4 2 (p) 
+ 15gt(p)) + 120D 3 (f!(p) - 23 ft(p)gl(p) + 23 f 2 4 (p)gt(p) - glip)) 
+9WD 2 (9fl(p)gl(p) - Mft(p)gi(p) + 9f 2 (p)gt(p))) 

-5D 3 f A (p)(3D - 16/ 2 (p))^ 4 (p)sin(6p 4 ))) (63) 

To simplify the expressions for the p4 action we define the expansion coefficients 
recursively and thus also list the odd expansion coefficients. However, after integra- 
tion over the momenta also in this case only even powers of \ijT contribute to the 
expansion of the pressure. Introducing further abbreviations, 

Sfj, = ^sin 2 (p„) and c k = -ic k , (64) 
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the expansion coefficients for the p4 action can be written as 



c = hip) (65) 

ci = -— ( - Si - S 2 - S 3 + 6S| + £>icos(2pi) + S 2 cos(2p 2 ) 

+S 3 cos(2p 3 ))sin(2p 4 ) (66) 

c 2 = tJtjO^ 2 + 6 ( 35 4 - Sisin 2 ( Pl ) - S 2 sin 2 (p 2 ) 

-S 3 sin 2 (p 3 )) sin 2 (p 4 ) - 2 cos 2 (p 4 ) (9S 4 

+ sin 2 (pi) ( - 3Si + sin 2 (p 4 )) + sin 2 (p 2 ) ( - 3S 2 + sin 2 (p 4 )) 

+ sin 2 (p 3 )(-3S 3 + sin 2 (p 4 )))) (67) 

c 3 = Ygf) ( 3j °ci (c 2 - 6c 2 ) - ( - 3 + cos(2pi) + cos(2p 2 ) 

+ cos(2p 3 )) cos 3 (p 4 ) sin(p 4 ) — 2 cos(p 4 ) sin(p 4 ) (12S 2 

+ sin 2 (pi) ( - 4Si + sin 2 (p 4 )) + sin 2 (p 2 ) ( - 4S 2 + sin 2 (p 4 )) 

+ sin 2 (p 3 )(-4S 3 + sin 2 (p 4 )))) (68) 

c 4 = ^^(6cos 4 (p 4 )(sin 2 (p 1 ) + sin 2 (p 2 ) + sin 2 (p 3 )) 

-3(3L>(c 4 - Ylc\c 2 + 12c 2 - 24cic 3 ) + 8( - 3S 2 + S x sin 2 ^) 

+S 2 sin 2 (p 2 ) + S 3 sin 2 (p 3 )) sin 2 (p 4 ) + ( - 3 + cos(2pi) 

+ cos(2p 2 ) + cos(2p 3 )) sin 4 (p 4 )) - 4 cos 2 (p 4 ) (l8S 2 

+ sin 2 (p!) ( - 6S1 + 11 sin 2 (p 4 )) + sin 2 (p 2 ) ( - 6S 2 + 11 sin 2 (p 2 )) 

+ sin 2 (p 3 ) ( - 6S 3 + 11 sin 2 (p 3 )))) (69) 

c 5 = - {3D (c? - 20c 3 c 2 - 60c 2 c 3 + 120c 2 c 3 + 60c 4 (c 2 + 2c 4 ) ) 

+20 ( — 3 + cos(2pi) + cos(2p 2 ) + cos(2p 3 )) cos 3 (p 4 ) sin(p 4 ) 

+8 cos(p 4 ) sin(p 4 ) (12S 2 + sin 2 (pi) ( - 4S X + 5 sin 2 (p 4 )) 

+ sin 2 (p 2 ) ( - 4S 2 + 5 sin 2 (p 4 )) + sin 2 (p 3 ) ( - 4S 3 + 5 sin 2 (p 4 )) ))(70) 

c 6 = — (6 - 9L>c? + 270L>c 4 c 2 - 1620£>c 2 c 2 + 1080.Dc 3 

6 6480L> V 1 1 12 2 

+1080 J Dc 3 c 3 - 6480Dc lC2 c 3 - 3240.Dc 2 - 3240D>c 2 c 4 
+6480D>c 2 c 4 - 6480Dcic 5 - 2 cos(2pi) - 2 cos(2p 2 ) - 2 cos(2p 3 ) 
+31 cos(2(p! - 2p 4 )) + 31 cos(2(p 2 - 2p 4 )) + 31 cos(2(p 3 - 2p 4 )) 
+24Si cos(2( Pl - p 4 )) + 24S 2 cos(2(p 2 - p 4 )) + 24S 3 cos(2(p 3 - p 4 )) 
-48Si cos(2p 4 ) - 48S 2 cos(2p 4 ) - 48S 3 cos(2p 4 ) + 288S^ cos(2p 4 ) 
-186 cos(4p 4 ) + 24Si cos(2(px + p 4 )) + 24S 2 cos(2(p 2 + p 4 )) 
+24S 3 cos(2(p 3 + Pi )) + 31 cos(2(p 1 + 2p 4 )) 

+31 cos(2(p 2 + 2p 4 )) + 31 cos(2(p 3 + 2p 4 ))). (71) 
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Note that we have defined here the coefficients q without iV r factors, which can be 
found in front of the integrals in (J53j) . 
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